using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// DQAWSE - Approximation to integral with algebraic and/or logarithmic singularities.
    /// </summary>
    [Serializable]
    public static class DqawseClass
    {
        /// <summary>
        /// Approximation to integral with algebraic and/or logarithmic singularities.
        /// </summary>
        /// <param name="f">Double precision function to be integrated.</param>
        /// <param name="a">Lower limit of integration.</param>
        /// <param name="b">Upper limit of integration.</param>
        /// <param name="alfa">The parameter in the weight function.</param>
        /// <param name="beta">The parameter in the weight function.</param>
        /// <param name="wgtfunc">Indicates which weight function is to be used: = 1:    (x-a)^alfa * (b-x)^beta; = 2:    (x-a)^alfa * (b-x)^beta * log(x-a); = 3:    (x-a)^alfa * (b-x)^beta * log(b-x); = 4:    (x-a)^alfa * (b-x)^beta * log(x-a) * log(b-x).</param>
        /// <param name="epsabs">The absolute accuracy requested.</param>
        /// <param name="epsrel">The relative accuracy requested.</param>
        /// <param name="abserr">The estimate of the modulus of the absolute error.</param>
        /// <param name="neval">The needed function evalations.</param>
        /// <param name="ier">The requested errorcode.</param>
        /// <returns>The approximation to the specified integral.</returns>
        public static double Dqawse(HardRealFunction f, double a, double b, double alfa, double beta, int wgtfunc,
                                    double epsabs, double epsrel, ref double abserr, ref int neval, ref int ier)
        {
            double[] alist = new double[QuadPackConstants.LIMIT];
            double[] blist = new double[QuadPackConstants.LIMIT];
            double[] rlist = new double[QuadPackConstants.LIMIT];
            double[] elist = new double[QuadPackConstants.LIMIT];
            double[] ri = new double[25];
            double[] rj = new double[25];
            double[] rh = new double[25];
            double[] rg = new double[25];
            double area;
            double area1;
            double area12 = 0;
            double area2;
            double a1;
            double a2;
            double b1;
            double b2;
            double centre;
            double errbnd;
            double errmax;
            double error1 = 0;
            double erro12;
            double error2 = 0;
            double errsum;
            double resas1 = 0;
            double resas2 = 0;
            double result;

            int[] iord = new int[QuadPackConstants.LIMIT];
            int iroff1;
            int iroff2;
            int k;
            int last;
            int limit;
            int maxerr;
            int nev = 0;
            ;
            int nrmax;

            limit = QuadPackConstants.LIMIT - 1;
            //  Test on validity of parameters. 
            ier = 6;
            neval = 0;
            rlist[0] = 0.0;
            elist[0] = 0.0;
            iord[0] = 0;
            result = 0.0;
            abserr = 0.0;
            if ((b <= a) || ((epsabs < 0.0) && (epsrel < 0.0)) || (alfa <= -1.0) || (beta <= -1.0) || (wgtfunc < 1) ||
                (wgtfunc > 4) || (limit < 1))
            {
                goto _999;
            }
            ier = 0;

            //  Compute the modified Chebyshev moments. 
            DqmomoClass.Dqmomo(alfa, beta, ri, rj, rg, rh, wgtfunc);

            //  Integrate over the invervals (a,(a+b)/2) and ((a+b)/2,b). 
            centre = 0.5 * (a + b);
            area1 = Dqc25sClass.Dqc25s(f, a, b, a, centre, alfa, beta, ri, rj, rg, rh, ref error1, ref resas1, wgtfunc,
                                       ref nev);
            neval = neval + nev;
            area2 = Dqc25sClass.Dqc25s(f, a, b, centre, b, alfa, beta, ri, rj, rg, rh, ref error2, ref resas2, wgtfunc,
                                       ref nev);
            neval = neval + nev;
            result = area1 + area2;
            abserr = error1 + error2;

            // Test on accuracy. 
            errbnd = (((epsabs) > (epsrel * Math.Abs(result))) ? (epsabs) : (epsrel * Math.Abs(result)));

            //  Initialization. 
            if (error1 >= error2)
            {
                alist[0] = a;
                alist[1] = centre;
                blist[0] = centre;
                blist[1] = b;
                rlist[0] = area1;
                rlist[1] = area2;
                elist[0] = error1;
                elist[1] = error2;
            }
            else
            {
                alist[0] = centre;
                alist[1] = a;
                blist[0] = b;
                blist[1] = centre;
                rlist[0] = area2;
                rlist[1] = area1;
                elist[0] = error2;
                elist[1] = error1;
            }
            iord[0] = 0;
            iord[1] = 1;
            if (limit == 1)
            {
                ier = 1;
            }
            if ((abserr <= errbnd) || (ier == 1))
            {
                goto _999;
            }
            errmax = elist[0];
            maxerr = 0;
            nrmax = 0;
            area = result;
            errsum = maxerr;
            iroff1 = 0;
            iroff2 = 0;

            //  Main loop. 
            for (last = 2; last < limit; last++)
            {
                a1 = alist[maxerr];
                b1 = 0.5 * (alist[maxerr] + blist[maxerr]);
                a2 = b1;
                b2 = blist[maxerr];

                area1 = Dqc25sClass.Dqc25s(f, a, b, a1, b1, alfa, beta, ri, rj, rg, rh, ref error1, ref resas1, wgtfunc,
                                           ref nev);
                neval = neval + nev;
                area2 = Dqc25sClass.Dqc25s(f, a, b, a2, b2, alfa, beta, ri, rj, rg, rh, ref error2, ref resas2, wgtfunc,
                                           ref nev);
                neval = neval + nev;

                //  Improve previous approximation and error test for accuracy. 
                area12 = area1 + area2;
                erro12 = error1 + error2;
                errsum += (erro12 - errmax);
                area += (area12 - rlist[maxerr]);
                if ((a == a1) || (b == b2))
                {
                    goto _30;
                }
                if ((resas1 == error1) || (resas2 == error2))
                {
                    goto _30;
                }

                //  Test for roundoff error. 
                if ((Math.Abs(rlist[maxerr] - area12) < (1.0e-5 * Math.Abs(area12))) && (erro12 >= (0.99 * errmax)))
                {
                    iroff1++;
                }
                if ((last > 9) && (erro12 > errmax))
                {
                    iroff2++;
                }
                _30:
                rlist[maxerr] = area1;
                rlist[last] = area2;

                //  Test on accuracy. 
                errbnd = (((epsabs) > (epsrel * Math.Abs(area))) ? (epsabs) : (epsrel * Math.Abs(area)));
                if (errsum <= errbnd)
                {
                    goto _35;
                }

                //  Set error flag in the case that number of intervals exceeds limit. 
                if (last == limit)
                {
                    ier = 1;
                }

                //  Set error flag in the case of roundoff error. 
                if ((iroff1 > 5) || (iroff2 > 19))
                {
                    ier = 2;
                }

                //  Set error flag in case of bad integrand behavior at interior points. 
                if ((((Math.Abs(a1)) > (Math.Abs(b2))) ? (Math.Abs(a1)) : (Math.Abs(b2))) <=
                    ((1.0 + 1.0e3 * Double.Epsilon) * (Math.Abs(a2) + 1.0e3 * Double.MinValue)))
                {
                    ier = 3;
                }
                //  Append the newly created intervals to the list. 
                _35:
                if (error2 <= error1)
                {
                    alist[last] = a2;
                    blist[maxerr] = b1;
                    blist[last] = b2;
                    elist[maxerr] = error1;
                    elist[last] = error2;
                }
                else
                {
                    alist[maxerr] = a2;
                    alist[last] = a1;
                    blist[last] = b1;
                    rlist[maxerr] = area2;
                    rlist[last] = area1;
                    elist[maxerr] = error2;
                    elist[last] = error1;
                }

                //  Call subroutine qsort to maintain the descending ordering. 
                DqsortClass.Dqsort(limit, last, ref maxerr, ref errmax, elist, iord, ref nrmax);

                //  Jump out of loop. 
                if ((ier != 0) || (errsum <= errbnd))
                {
                    break;
                }
            }
            result = 0.0;
            for (k = 0; k <= last; k++)
            {
                result += rlist[k];
            }
            abserr = errsum;
            _999:
            return result;
        }
    }
}